Mock min and max temp list for testing

t_min_list = list("1"=-12, "2"=-8, "3"=-4, "4"=-3, "5"=8, "6"=10, "7"=11, "8"=15, "9"=18, "10"=3, "11"=-2, "12"=-10)

t_max_list = list("1"=-1, "2"=-3, "3"=-2, "4"=0, "5"=18, "6"=22, "7"=28, "8"=23, "9"=22, "10"=6, "11"=3, "12"=-7)

library(ggplot2)
library(viridis)

## Loading required package: viridisLite

library(ggExtra)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

PLOT

createTimeSeries <- function (data, y_lab) {

  # Most  plot
  p <- ggplot(data, aes(x=Year, y=Value, group = 1)) +
    geom_line(color="steelblue") +
    xlab("Year") +
    ylab(y_lab)
  
#  +
#    theme_ipsum() +
#    theme(axis.text.x=element_text(angle=60, hjust=1))
  
  p <- ggplotly(p)
  
  p
}

NFFD

nffd_param <- read.csv(file = "../optimizedParameterTables/param_NFFD.csv", sep=',', header = TRUE)

# nffd: number of frost free days
# m: month of the year
# tm: min temperature for that month
nffd <- function(m, tm) {
  
  optimized_params <- nffd_param[nffd_param$Month == m,]

  a <- optimized_params$a
  b <- optimized_params$b
  t0 <- optimized_params$T0
  
  cat(t0)

  return( a/(1 + exp(-(tm - t0)/b)))
}

nffd(2, -3)

## 3.8

## [1] 0.03910693

Fucntion to compute NFFD with data

compute_nffd <- function (month, scenario, model) {

  month_two_digits = sprintf("%02d",month)
  file = paste("../data/ensmax.COM.Tmin", month_two_digits, ".csv", sep="")
  
  tmin_month <- read.csv(file, sep=',', header = TRUE)

  tmin_month <- tmin_month[tmin_month$scenario == 'historical',][c("Year",model)]
  
  nffd_value <- lapply(tmin_month[model], function(tm) nffd(month, tm))
  nffd_year <- tmin_month$Year
  #nffd_date <- as.Date(with(tmin_month,paste(Year, month, "01", sep="-")),"%Y-%m-%d")
  
  nffd_df <- data.frame(nffd_year,nffd_value)
  names(nffd_df) <- c("Year", "Value")
  
  return(nffd_df)
}

Test Monthly NFFD

scenario = "historical"
model =  "ACCESS.ESM1.5"

jan_nffd <- compute_nffd(1, scenario, model)

## 3.82

feb_nffd <- compute_nffd(2, scenario, model)

## 3.8

mar_nffd <- compute_nffd(3, scenario, model)

## 3.6

apr_nffd <- compute_nffd(4, scenario, model)

## 3.13

may_nffd <- compute_nffd(5, scenario, model)

## 2.78

jun_nffd <- compute_nffd(6, scenario, model)

## 2.24

jul_nffd <- compute_nffd(7, scenario, model)

## 1.23

aug_nffd <- compute_nffd(8, scenario, model)

## 1.71

sep_nffd <- compute_nffd(9, scenario, model)

## 2.72

oct_nffd <- compute_nffd(10, scenario, model)

## 3.23

nov_nffd <- compute_nffd(11, scenario, model)

## 3.47

dec_nffd <- compute_nffd(12, scenario, model)

## 3.63

createTimeSeries(jan_nffd, "NFFD Historical Jan")

createTimeSeries(feb_nffd, "NFFD Historical Feb")

createTimeSeries(mar_nffd, "NFFD Historical Mar")

createTimeSeries(apr_nffd, "NFFD Historical Apr")

createTimeSeries(may_nffd, "NFFD Historical May")

createTimeSeries(jun_nffd, "NFFD Historical Jun")

createTimeSeries(jul_nffd, "NFFD Historical Jul")

createTimeSeries(aug_nffd, "NFFD Historical Aug")

createTimeSeries(sep_nffd, "NFFD Historical Sep")

createTimeSeries(oct_nffd, "NFFD Historical Oct")

createTimeSeries(nov_nffd, "NFFD Historical Nov")

createTimeSeries(dec_nffd, "NFFD Historical Dec")

Function to compute seasonal NFFD

#Winter
season_nffd <- function(months, scenario, model) {
  
  list_of_monthly_nffd = list()
  
  for(month in months) {
      list_of_monthly_nffd[[month]] <- compute_nffd(month, scenario, model) 
  }

  
  nffd_all <- data.table::rbindlist(list_of_monthly_nffd)
  agg_summ_nffd_all = aggregate(Value~Year,nffd_all,FUN = sum)

  return(agg_summ_nffd_all)
}

Compute seasonal NFFD’s

scenario = "historical"
model =  "ACCESS.ESM1.5"

winter_months = c(12,1,2)
nffd_winter <- season_nffd(winter_months,scenario, model)

## 3.633.823.8

spring_months = c(3,4,5)
nffd_spring <- season_nffd(spring_months,scenario, model)

## 3.63.132.78

summer_months = c(6,7,8)
nffd_summer <- season_nffd(summer_months,scenario, model)

## 2.241.231.71

autumn_months = c(9,10,11)
nffd_autum <- season_nffd(autumn_months,scenario, model)

## 2.723.233.47

createTimeSeries(nffd_winter, "NFFD Historical Winter")

createTimeSeries(nffd_spring, "NFFD Historical Spring")

createTimeSeries(nffd_summer, "NFFD Historical Summer")

createTimeSeries(nffd_autum, "NFFD Historical Autum")

FFP, bFFP and eFFP

# bffp: the day of the year on which FFP begins
# t_min_list: named list of monthly minimum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
# nffd: number of frost-free days.
bffp <- function(td, nffd, t_min_list) {
  
  tmin4 <- t_min_list[["4"]]
  tmin6 <- t_min_list[["6"]]
    
  return(352.1358994 + -0.021715653 * tmin4^2 + -3.542187618 * tmin6 + 0.020359471 * tmin6^2 - 4.897998097 * td + 0.033521327 * td^2 - 2.164862277 * nffd + 0.006767633 * nffd^2 - 0.00000929 * nffd^3 + 0.043516586 * (td * nffd) - 0.00000253 * (td * nffd)^2)
}

# effp: the day of the year on which FFP ends 
# t_min_list: named list of monthly minimum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
effp <- function(nffd, t_min_list) {
  
  tmin9 <- t_min_list[["9"]]
  tmin10 <- t_min_list[["10"]]
  tmin11 <- t_min_list[["11"]]
    
  return(243.7752209 + 4.134210825 * tmin9 - 0.162876448 * tmin9^2 + 1.248649021 * tmin10 + 0.145073612 * tmin10^2 + 0.004319892 * tmin10 + -0.005753127 * tmin10^2 - 0.06296471 * nffd + 0.000399177 * nffd^2)
  
}

ffp <-function(effp,bffp) {
  return(effp - bffp)
}

Test

bffp_test <- bffp(10, jun_nffd, t_min_list)
effp_test <- effp(jun_nffd, t_min_list)
ffp_test <- ffp(effp_test, bffp_test)

PAS

pas_param <- read.csv(file = "../optimizedParameterTables/param_PAS.csv", sep=',', header = TRUE)

# pas: precipitation as snow
# m: month of the year
# tm: min temperature for that month
pas <- function(m, tm) {
  
  optimized_params <- pas_param[pas_param$Month == m,]

  b <- optimized_params$b
  t0 <- optimized_params$T0

  return( 1/(1 + exp(-(tm - t0)/b)))
}

Test

pas(5, 10)

## [1] 0.9998033

EMT, EXT

# emt: extreme minimum temperature
# t_min_list: named list of monthly minimum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
emt <- function(t_min_list, td) {
  
  tmin1 <- t_min_list[["1"]]
  tmin12 <- t_min_list[["12"]]

  # tminx: minimum temperature over the year
  tminx <- min(sapply(t_min_list, min))

  
  return(-23.02164 + 0.77908 * tmin1 + 0.67048 * tmin12 + 0.01075 * tminx^2 + 0.11565 * td)
}

# ext: extreme maximum temperature
# t_max_list: named list of monthly maximum temperature for each month
# td: difference between the mean warmest monthly temperature and the mean coldest monthly temperature
ext <- function(t_max_list, td) {
  
  tmax7 <- t_max_list[["7"]]
  tmax8 <- t_max_list[["8"]]

  # tmaxx: maximum temperature over the year
  tmaxx <- max(sapply(t_max_list, max))

  
  return(10.64245 + -1.92005 * tmax7 + 0.04816 * tmax7^2 + 2.51176 * tmax8 - 0.03088 * tmax8^2 - 0.01311 * tmaxx^2 + 0.33167 * td - 0.001 * td^2)
}

Test

emt(t_min_list, 23)

## [1] -34.86745

RH

# es: saturated vapour pressure at a temperature t
# t: air temperature
es <- function(t) {
  
  svp <- 0.6105 * exp((17.273*t)/(t+237.3))
  
  if(t < 0) {
    return(svp*(1 + (t*0.01)))
  } else {
    return(svp)
  }
}

# *******************************************
# Question, should this be tmin or tmin_mean?
# *******************************************

# rh: relative humidity
# tmin_mean: monthly mean minimum air temperature
# tmax_mean: monthly mean maximum air temperature
rh <- function(tmin_mean, tmax_mean) {
  es_avg = (es(tmin_mean)+ es(tmax_mean))/2
  
  return((100 * es(tmin_mean)/es_avg))
}

es(-10)

## [1] 0.2569797

rh(-12, 3)

## [1] 44.07521

DD

dd_param_below_0 <- read.csv(file = "../optimizedParameterTables/param_DD_S1.csv", sep=',', header = TRUE)
dd_param_above_5 <- read.csv(file = "../optimizedParameterTables/param_DD_S2.csv", sep=',', header = TRUE)
dd_param_below_18 <- read.csv(file = "../optimizedParameterTables/param_DD_S3.csv", sep=',', header = TRUE)
dd_param_above_18 <- read.csv(file = "../optimizedParameterTables/param_DD_S4.csv", sep=',', header = TRUE)

dd <- function(m, tm) {
  
# ***********************************************
# Question, I think it should be dd < 5 (not > 5)
# ***********************************************
  dd_param <- ''

  if(tm < 0) {
    dd_param <- dd_param_below_0
  } else if(tm < 5) {
    dd_param <- dd_param_above_5[dd_param_above_5$Region == "All"]
  } else if(tm < 18) {
    dd_param <- dd_param_below_18
  } else {
    dd_param <- dd_param_above_18
  }
  
  optimized_params <- dd_param[dd_param$Month == m,]
    
  k <- optimized_params$a
  a <- optimized_params$a
  b <- optimized_params$b
  t0 <- optimized_params$T0
  c <- optimized_params$a
  beta <- optimized_params$a
  
  if(tm > k) {
    return( a/(1 + exp(-(tm - t0)/b)))
  } else {
    return(c + (beta * tm))
  }
}

dd(4, 14)

## [1] 4878.885

Test functions